Method for determining the accuracy of measured environmental data parameters

ABSTRACT

A method of estimating expected errors of environmental data parameters based on radiance measurements obtained from visible infra red radiometric satellite sensors (VIIRS) orbiting the earth, comprising the steps of: obtaining N radiance measurements of a surface body I I , . . . ,I N  defining matrix I depending on p unknown surface and atmospheric parameters T i , . . . T p , defining a matrix T; generating a forward model I=f(T) for obtaining the I radiance measurements from the p parameters; choosing an initial set of values for the p parameters and linearizing f(T) about the initial values to obtain a linearized forward I=s+Hθ as f i  (T 01 , T 02 , . . . T 0p )+ΣH ij  θ j  where I=1,2, . . . N and θ j  =T oj  and H ij  =δf i  /δT j  and where θ is a column matrix θ l , . . . θ p  and H is a matrix of H ij  values; adding measurement noise vector w of noise values to the forward model; determining the covariance of the measurement noise w to obtain a covariance matrix C; and manipulating the matrices H and C according to the equation C EDR  =(H T  C -1  H) -1  to obtain matrix element C EDR  indicative of the expected errors in the values of T 1 , . . . T p  parameters.

FIELD OF THE INVENTION

The present invention relates generally to satellite radiometry measurements for measuring environmental data parameters, and more particularly, to a method for determining the accuracy with which instruments providing said radiometry measurements yield the true values of the measured parameters at the earth's surface.

BACKGROUND OF THE INVENTION

There are certain phenomena on our earth which have great significance from military and/or commercial points of view. One such phenomenon is the water temperature of the oceans or other large bodies of water. Knowing the temperature of the water over large areas is of importance for military purposes, as well as from a commercial point of view, since ocean water temperature affects commercial fishing. Theoretically, one can take direct measurements of the water temperature at specific locations. However, due to the large surface area of the oceans, taking direct water temperature measurements over any large surface area of hundreds of thousands of square miles is clearly impractical.

Similarly, cloud measurements in general, and more particularly, ice particle size measurements of cirrus clouds, are also important in imaging certain bodies in order to remove or limit interfering clouds from any satellite images, as well as to obtain measurements of precipitation occurring within the atmosphere. Satellite sensor instruments have been used to measure specified radiance of certain ground-based bodies to within a specified error. That is, the satellite sensors themselves make known radiance measurements at the top of the atmosphere of bodies located at the earth's surface, such as water temperature.

However, the goal is to determine a given Environmental Data Record (EDR) or surface parameter to within a specified error. It is therefore crucial to relate the radiance measurement errors, which we know how to predict, to the expected errors in the EDR's. For many EDR's it is relatively easy to create forward models which, given specific values of the EDR parameters, predict the radiance leaving the top of the atmosphere. What is not so easy, however, is to exercise these models over a range of environmental parameters large enough to find the level of radiance error which produces adequate EDR estimates. Fortunately, as will be shown herein, signal processing formulas applied to linearized versions of an EDR forward model, permit one to proceed directly from the known radiance error (i.e. the satellite sensor noise of the point design) to the expected errors in the EDR. These formulas also specify algorithms, based on the linearized forward model, which can be used to estimate EDR's from their associated radiances.

As will be shown in the detailed description below, a basic statistical methodology for connecting EDR errors to radiance measurement errors will be presented in Section 1, followed by two examples in Sections 2 and 3 showing how to deal with VIIRS-based EDR's when the forward models can be written down algebraically. The Sea-Surface Temperature EDR analyzed in Section 2 shows that the predicted EDR error using the method disclosed herein matches NOAA's past experience. This provides support that the statistical methodology works as expected. The cirrus cloud EDR analyzed in Section 3 shows how the methodology can be used to analyze the performance different versions of an EDR algorithm. Section 4 shows how to handle EDR's where only a computerized version of the forward model (such as MODTRAN or FASCODE) exists. The technique described in Section 4 applies equally well, of course, to those VIIRS-based EDR's where only computerized versions of the forward model are easily available.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide a method for estimating expected errors of environmental data parameters based on radiance measurements obtained from visible infra red radiometric satellite sensors (VIIRS) orbiting the earth, comprising the steps of obtaining N radiance measurements of a surface body I_(I), . . . I_(N) defining matrix I depending on p unknown surface and atmospheric parameters T_(i), . . . T_(p), defining a matrix T; generating a forward model I=f(T) for obtaining the I radiance measurements from the p parameters, where f(T) comprises a matrix; choosing an initial set of values for the p parameters and linearizing f(T) about the initial values to obtain a linearized forward I=s+Hθ as f_(i) (T₀₁, T₀₂, . . . T_(0p))+ΣH_(ij) θ_(j) where I=1,2, . . . N and θ_(j) =T_(j) -T_(oj) and H_(ij) =δf_(i) /δT_(j) and where θ is a column matrix θ₁, . . . θ_(p) and H is a matrix of H_(ij) values; adding measurement noise vector w of noise values to the forward model; determining the covariance of the measurement noise w to obtain a covariance matrix C; and manipulating the matrices H and C according to the equation C_(EDR) =(H^(T) C⁻¹ H)⁻¹ to obtain matrix element C_(EDR) indicative of the expected errors in the values of T₁, . . . T_(p) parameters.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a standard two-layer transmission-emission model for the ocean surface and atmosphere according to the present invention.

FIG. 2 illustrates the transmission curves as a function of wavenumber listed in Table 1 and shows how well they match the straight-line approximations generated by the values of M and A_(atm).

FIG. 3 represents a graph of functions P and α for T=290 deg K., 300 deg K., 310 deg K. and 800 cm⁻¹ ≦v≦1000 cm⁻¹.

FIG. 4 shows the relative error involved in approximating P by α for three different temperatures.

FIG. 5 illustrates the physical model for measuring cirrus cloud temperature and ice-particle size using radiance measurements according to the present invention.

FIG. 6 provides an overall flow chart for estimating the expected errors of environmental data parameters based on radiance measurements from VIIRS sensors according to the present invention.

FIG. 7 is a flow chart depicting the steps for predicting errors in sea surface temperature measurements according to the present invention.

FIGS. 8A-B represent a flow chart depicting a method for estimating error involved in determining ice particle size of cirrus clouds using VIIRS radiance measurements according to the present invention.

FIG. 9 contains Table 1 which provides a set a of ratios for midlatitude winter to midlatitude summer and for tropical to midlatitude summer temperature values a function of wavenumber.

FIG. 10 contains Table 2 which provides a ration of the slope to the y-intercept from the value of Table 1.

FIG. 11 contains Table 3 which lists the expected errors in sea-surface temperature ΔT, based on the VIIRS point design, for different values of the sea-surface temperature.

FIG. 12 contains Table 4 which provides a list of expected sea-surface errors based on AVHRR specifications for different values of sea-surface temperature.

FIG. 13 contains Table 5 which illustrates parameter values representing bands 3 and 4 of the VIIRS point design for the 3.7 μm bands of the 10-12 μm bands.

FIG. 14 contains Table 6 which represents ice particle size measurement values, and the absolute and relative errors associated with those measurements over a range of T_(C) temperature values.

DETAILED DESCRIPTION OF THE INVENTION

Section 1: The Statistical Methodology

With reference now to FIG. 6, as part of the statistical methodology for VIIRS-based EDR parameters, we suppose there exist N measurements (e.g. N in-band sensor signals) which depend on p unknown atmospheric and surface parameters. In order to make the mathematical notation compact the N measurements I₁, I₂, . . . , I_(N) are stacked in a column and treated as a column vector I. ##EQU1##

The p unknown surface and atmospheric parameters T₁, T₂, . . . , T_(p) are treated similarly. ##EQU2##

In general the T are the EDR parameters which are to be determined from the satellite sensor measurements I. For example, there may be N=2, p=2 for a sea-surface temperature EDR, where I₁,2 are two radiance measurements in the 10-12 μm atmospheric window, and T₁,2 are respectively the sea-surface temperature and an atmospheric transmission-emission parameter as shown in equations (2.8a-d) below. Note that for CrIS-based atmospheric profiles there may be N≅1500 for the spectral channels and p≅70 for a 35-layer atmospheric model with temperature and water-vapor concentrations to be determined for each of the layers.

The forward model going from T to I can be written as

    I=f(T)

where f(T) stands for a column vector of functions all depending on T. ##EQU3## The next step is to choose a standard set of values for the surface and atmospheric parameters. ##EQU4## and linearize f₁, f₂, . . . , f_(N) about these values. ##EQU5## where

    i=1,2, . . . , N,

    θ.sub.j =T.sub.j -T.sub.0j,                          (1.3c)

and ##EQU6## If one defines ##EQU7## and

    H=matrix of H.sub.ij values                                (1.4c)

then the linearized forward model in equation (1.2a) can be written as

    I=Hθ+s                                               (1.5)

Measurement noise associated with the N radiance measurements is represented as a vector w of N random numbers w₁, W₂, . . . , W_(N) having a jointly gaussian probability distribution.

The probability that w₁ has a value between w₁ and w₁ +dw₁, w₂ has a value between w₂ and w₂ +dw₂, . . . , w_(N) and W_(N) +dw_(N) is

    p(w.sub.1, w.sub.2, . . . w.sub.N)dw.sub.1 dw.sub.2. . . dw.sub.N

where ##EQU8## Note that this is a zero-mean gaussian distribution. In equation (1.7a) C is the covariance matrix of the w₁, W₂, . . . , w_(N) random numbers, det C is the determinant of C, and w^(T) =(w₂ w,._(N) w ) is the transpose of w (which is, of course, a row vector). The ij'th component of matrix C is

    (C).sub.i =C.sub.i =E((w.sub.i -E(w.sub.i))·(w.sub.j -E(w.sub.j)))=E(w.sub.i w.sub.j)                          (1.7b)

where E stands for the usual expectation operator. In the last step of equation (1.7b) we use that E(w_(i))=0 for all values of I because the w_(i) obey a zero-mean gaussian distribution. For any two random numbers r₁, r₂ (which do not necessarily obey a zero mean distribution) we have

    E(r.sub.1,2)=mean value of r.sub.1,2                       (1.7c)

    E((r.sub.1,2 -E(r.sub.1,2)).sup.2)=variance of r.sub.1,2   (1.7d)

    E((r.sub.1 -E(r.sub.1)).(r.sub.2 -E(r.sub.2)))=covariance of r.sub.1 and r.sub.2                                                   (1.7e)

One can then construct a complete linearized forward model, including measurement noise, by adding w, the vector or random measurement errors, to the right-hand side of equation (1.5) to obtain:

    I=Hθ+s+wtm (1.8)

The error in the estimate of θ₁, θ₂, . . . , θ_(p) is then given by the covariance matrix

    C.sub.EDR =(H.sup.T I.sup.-1 H).sup.-1                     (1.9a)

where the "T" superscript denotes the matrix transpose and the "-1" superscript denotes the matrix inverse. The covariance matrix C_(EDR) follows from the theorem established in Steven M. Kay's, Fundamentals of Statistical Signal Processing: Estimation Theory, page 97, PTR Prentiss-Hall, Inc., Englewood Cliff, N.J., 1993 and incorporated herein by reference. The matrix elements of C_(EDR) specify the expected errors in the values of θ₁,θ₂, . . . θ_(p), which is the same as the expected errors in the values of T₁,T₂, . . . T_(p) because the θ's are just the differences of the p surface and atmospheric parameters from their standard values (see equation (1.3c)). The minimum-variance, unbiased estimator of θ given noise-contaminated values of I0 (the noise amplitude is specified by the covariance matrix C) and known values of s is

    θ.sub.estimate =(H.sup.T C.sup.-1 H).sup.-1 H.sup.T C.sup.-1 (I-s)(1.91b)

Equation (1.9a) specifies the variance, that is, the error, in the estimate of θ given by equation (1.9b). If the forward model in equation (1.8) is strictly true rather than a linearized version of equation (1.2a), then equation (1.9b) is the best possible way to estimate θ and equation (1.9a) gives the lowest possible error with which θ can be estimated.

The constant s in equation (1.8) and (1.9b) is a known quantity so we can always define a new set of measurements D such that

    D-I-s                                                      (1.10a)

and

    D=Hθ+w                                               (1.10b)

If θ has been measured often enough in the past that one can estimate the probability that the next value of θ will take on any particular value, and if it is known that this probability distribution for θ looks gaussian with a covariance matrix C.sub.θ such that

    (C.sub.0).sub.ij =E((θ.sub.i -E(θ.sub.i))(θ.sub.j -E(θ.sub.j))).                                      (1.11a)

then it is also known that the expected error in our estimate of θ will now be

    C.sub.EDR =(C.sub.θ.sup.-1 +H.sup.T C.sup.-1 H).sup.-1(1.11b)

when using the estimation formula (see page 391 of Fundamentals of Statistical Processing: Estimation Theory)

    θ.sub.estimate =E(θ)+[(C.sub.0.sup.-1 +H.sup.T C.sup.-1 H).sup.-1 H.sup.T C.sup.-1 ].(D-H.E( θ))            (1.11c)

If the D values are unknown so as to be unable to provide information regarding the values of θ, one may proceed by obtaining a best estimate as to the true values of θ as E(θ)=(E(θ₁), E(θ₂), . . . , E(θ₃))^(T), the mean of the set of previously-determined θ values. However, we do in fact have the θ values generated by the true θ values, the right-hand side of equation (1.11c) uses the information provided by the D values to adjust the already known values of E(θ) to get a better estimate of θ. In equations (1.3a,b,c) above, we can choose to linearize the forward model about a T₀ vector such that E(θ)=0 (indeed that would be the most natural way of linearizing the forward model). This simplifies equation (1.11c) to

    θ.sub.estimate =[(C.sub.θ.sup.-1 +H.sup.T C.sup.-1 H) .sup.-1 H.sup.T C.sup.-1 ]. D                                     (1.11d)

Equation (1.11 d) represents a standard linear algorithm used to retrieve temperature and water vapor from a set of FTIR spectra. The algorithm of this type is disclosed in pages 12-14 of J. P. Kerekes, A Retrieval Error Analysis Technique for Passive Infrared Atmospheric Sounders, Lincoln Laboratory MIT, Technical Report 978, Jul. 8, 1993, ESC-TR-92-204, incorporated herein by reference. Note that when (1.11d) is used, equation (1.11b) automatically gives us the expected error in the retrieved temperature and water-vapor profiles. If no information is available about which temperatures and water-vapor concentrations are most likely to occur, then all the elements of C.sub.θ are zero and equations (1.11b,d) reduce to (1.9a,b).

Section 2: Sea-Surface Temperature

A forward model for determining the sea-surface temperature from satellite observations of a cloud-free region of the ocean is given as

    I(v,atm)=B(v,T.sub.s)·τ(v,atm)+B(v,T.sub.a)·[1-τ(v,atm)]                                                    (2. 1)

where I is the radiance leaving the earth's atmosphere (in watt/cm² /sr/cm⁻¹), T_(s) is the sea-surface temperature (in deg K), T_(a) is an effective temperature (in deg K) describing the atmospheric self-radiance in the 10-12 μm window region, τ is the dimensionless atmospheric transmission from the ocean surface to the top of the atmosphere, v is a wavenumber (in cm⁻¹) in the 10-12 μm window region, "atm" is shorthand for a long list of atmospheric parameters determining the value of the transmittance τ at wavenumber v, and B(v,T) is the Planck function

    B(v,T)=C.sub.1 v.sup.3 P(C.sub.2 vT.sup.-1)                (2.2a)

where ##EQU9##

Equation (2.1) represents a standard two-layer transmission-emission model for the ocean surface and atmosphere as shown in FIG. 1. Empirical data has shown that T_(a) varies by less than 1 deg K in the 10-12 μm window region, so the dependence of T_(a) on v can be neglected in equation (2.1). Another good approximation is that in the 10-12 μm window

    τ(v,atm)≅1-A.sub.atm ·K(v)          (2.3)

where K is a function only of wavenumber and the effect of all the atmospheric parameters can be lumped into a single constant A_(atm). As a check to equation (2.3), a MODTRAN simulation run for 35 wavenumbers between 820 cm⁻¹ and 990 cm⁻¹ has been performed using three different standard atmospheres: tropical, midlatitude summer, and midlatitude winter. Table 1 (FIG. 9) gives the ratios of [1-τ(v,atm)] for midlatitude winter to midlatitude summer and for tropical to midlatitude summer, and we see that while some correlation exists with the approximation given in (2.3), there are definitely wavenumbers where it does not work very well. We will, however, be interested in wide-bandwidth radiance signals, so we should also examine how well equation (2.3) works as a moving average relationship. Performing a linear regression of [1-τ(v,atm)] against v for the values of 1-τ and v given in Table 1 yields the ratio of the slope to the y-intercept as shown in Table 2 (FIG. 10). As one can ascertain, the ratios are very similar for all three atmospheres (in fact they are equal to within 4%). Section 5 shows that this is exactly the result to be expected if equation (2.3) held true, so equation (2.3) can be used with some confidence as long as it is applied only to wide-bandwidth radiance signals. FIG. 2 gives the transmission curves τ(v,atm) listed in Table 1 and shows how well they match the straight-line approximations generated by the values of M and A_(atm) in equations (2.4a,b).

    τ(v,atm)≅1-A.sub.atm ·K.sub.L (v)   (2.4a)

where

    K.sub.L (v)=1-M·v                                 (2.4b)

for

    M=8.0×10.sup.-4 cm,

1.8 for standard tropical atm

and

A_(atm) =1.2 for standard midlatitude summer atm

0.27 for standard midlatitude winter atm

Approximating the function K(v) in equation (2.3) by the straight line K_(L) (v) given in (2.4b) and substitute (2.4a,b) into equations (2.1) yields.

    I(v,atm)=B(v,T.sub.s)+K.sub.L (v)·A.sub.atm ·[B(v,T.sub.a)-B(v,T.sub.s)]                     (2.5)

From page 5042 of C. Prabhakara, G. Dalu, and V. G. Kunde's article on "Estimation of Sea Surface Temperature From Remote Sensing in the 11- to 13 μm Window Region", J. Of Geophysical Research, vol. 79, No. 33, 1974, pp. 5039-5044, incorporated herein by reference, we know that T_(s) and T_(a) are almost always within 10 deg K of each other, so for a given patch of ocean at a specified time of year, the Planck function B(v,T) will only need to be evaluated over a limited range of V and T values. A reasonable range for midlatitude summer is

    285 deg K≦T≦315 deg K                        (2.6a)

    800 cm.sup.-1 ≦v≦1000 cm.sup.-1              (2.6b)

For this range of v and T function P(C₂ vT⁻¹) can be approximated as

    P(C.sub.2 vT.sup.-1)≅α(C.sub.2 vT.sup.-1)=Ω·(C.sub.2 vT.sup.-1).sup.-n    (2.7a)

where

    Ω=8.2527                                             (2.7b)

and

    n=4.4                                                      (2.7c)

FIG. 3 shows functions P and α for T=290 deg K, 300 deg K, 310 deg K and 800 cm⁻¹ ≦v≦1000 cm⁻¹ ; one can see that a does a reasonable job of approximating P. FIG. 4 shows the relative error involved in approximating P by α, and for the three temperatures graphed in FIG. 3, the relative error is 3% or less. We substitute (2.2a) and (2.7a) into equation (2.5) to get

    I(v,atm)=C.sub.b v.sup.-3-n θ.sub.s +C.sub.B K.sub.L (v)v.sup.3-n θ.sub.a                                             (2.8a)

where

    C.sub.B =ΩC.sub.1 C.sub.2.sup.-n ≅1.98×10.sup.-12 watt/cm.sup.2 /sr/(cm.sup.-1).sup.4 /(cm-degK).sup.n,     (2.8b)

    θ.sub.s =T.sub.s.sup.n                               (2.8c)

and

    θ.sub.a =A.sub.atm ·(T.sub.a.sup.n -T.sub.s.sup.n)(2.8d)

Equations (2.8a-d) reduce the nonlinear problem of equation (2.5) into a linear problem in two unknowns: θ_(s) the n'th power of the sea-surface temperature, and θ, an atmospheric variable describing the combined effect of T_(s), T_(a), and A_(atm). Therefore, for this particular VIIRS-based EDR we do not have to linearize the forward model to use equations (1.9a,b), although that could have been done had it turned out to be necessary.

The VIIRS point design specifies two bandwidths in the 10-12 μm atmospheric window which are labeled band 1 and band 2. Band 1 ranges from wavenumbers V_(1a) to v_(1b) and band 2 ranges from v_(2a) to v_(2b) where

    v.sub.1a =800 cm.sup.-1, v.sub.1b =870 cm.sup.-1           (2.9a)

    V.sub.2a =885 cm.sup.-1, v.sub.2b =971 cm.sup.-1           (2.9b)

The noise amplitudes in band 1 and 2, based on the point design specification of NedT=0.1 deg K at T=293 deg K, are σ₁ and σ₂ respectively where

    σ.sub.1 =1.211×10.sup.-6 watt/cm.sup.2 /sr     (2.10a)

    σ.sub.2 =1.284×10.sup.-6 watt/cm.sup.2 /sr     (2.10b)

One obtains σ₁ and σ₂ from the formula (I=1,2) ##EQU10## where the second step in equation (2.10c) comes from the approximation for B(v,T) given by equations (2.2a,b) and (2.7a-c).

Integrating equation (2.8a) over bands 1 and 2 gives (I=1,2) ##EQU11## where equation (2.4b) has been used to go from (2.11a) to (2.11b). Equation (2.1 lb) is written in matrix form as

    I=Hθ                                                 (2.12a)

where ##EQU12## The forward model with additive noise w=(w₁, w₂)^(T), where w₁ is the random error in band 1 and w₂ is the random error in band 2, can be written as

    I=Hθ+w                                               (2.13a)

From equations (2.10a,b) we see that, assuming gaussian signal noise with independent noise pulses in bands 1 and 2, the noise covariance matrix representing the measurement errors specified by w is ##EQU13## Equation (2.13a) is identical to equation (1.8) in Section 1 above (with s=0). Therefore, the C_(EDR) covariance matrix may be written as ##EQU14## The part of C_(EDR) of interest to the present application is the "1,1" component specifying the expected variance in the sea-surface variable θ_(s) =T_(s) ^(n). The error in T_(s), called .increment.t_(s), comes from ##EQU15## Equations (2.7c) and (2.15 b) are used in Table 3 (FIG. 11) to list .increment.T, for different values of the sea-surface temperature. The values in Table 3 run from 0.47 deg K to 0.59 deg K and are a good match to NOAA's past experience with AVHRR data (typical errors are 0.5 deg K to 0.7 deg K, see reference 14). This is not surprising because the VIIRS point design closely resembles the AVHRR specifications in the 10-12 μm atmospheric window. Repeating our calculations for the AVHRR bands in the 10-12 μm atmospheric window using AVHRR noise specifications,

    σ.sub.1,AVHRR =1.470×10.sup.-6 watt/cm.sup.2 /sr

    σ.sub.2,AVHRR =1.740×10.sup.-6 watt/cm.sup.2 /sr

the results show that the expected temperature errors range from 0.61 deg K to 0.78 deg K (see Table 4, FIG. 12). Again, there exists a reasonable match to NOAA's past experience in retrieving sea-surface temperatures. As shown in section 6, an algorithm may be implemented to calculate the temperature errors as shown in Tables 3 and 4. A flow chart depicting the above operation is illustrated in FIG. 7.

Section 3: Ice Particle Size in Cirrus Clouds

Here there is shown that the statistical methodology can be used to investigate the accuracy of different versions of an EDR algorithm. In this analysis, the processing steps involve directly linearizing the forward model to obtain an accurate error estimation rather than re-arranging into linear format as was done in the Section 2 example above. Note that the steps involved in obtaining an accurate error estimate may be implemented by means of a computer program resident on a computer host such as SPARC station or other commercially available computer systems. A flow chart depicting this method is illustrated in FIG. 8.

As is well known, the phenomenological relation between cirrus cloud temperature and ice-particle size is ##EQU16## where

L=average ice-particle size (in μm)

T_(c) =cirrus=cloud temperature (in deg K)

c₀ =326.3 μm

c₁ =12.42 μm/deg K

c₂ =0.197 μm/(deg K)²

c₃ =0.0012 μm/(deg K)³

Equation (3.1a) determines the ice-particle size once the cirrus clouds temperature is known. The forward model for the cirrus cloud temperature T_(c) is

    I.sub.3 (v)=(1-ε.sub.3 (v))·I.sub.G3 (v)+ε.sub.3 (V)·B(v,T.sub.c)                                 (3.2a)

for v a wavenumber (in cm⁻¹) in the 3.7 μm region of the earthshine radiance spectrum, I₃ the radiance (in watt/sm² /sr/cm)¹ in the 3.7 μm region leaving the cloud, IG_(G3) the radiance (in watt/cm² /sr/cm⁻¹) in the 3.7 μm region coming up from the earth beneath the cloud, ε₃, the dimensionless cloud emissivity in the 3.7 μm region, and B the Planck function defined above in equations (2.2a,b). A similar equation can be written for the radiance in the 10-12 μm region of the earthshine spectrum.

    I.sub.4 (v)=(1-ε.sub.4))·I.sub.G4 (v)+ε.sub.4 (v)·B(v,T.sub.c)                                 (3.2b)

with I₄, I_(G4), and ε₄ the same as I₃, I_(G3), and E₃ except that they apply to the 10-12 μm instead of the 3.7 μm spectral region. FIG. 5 shows that the conceptual basis of equations (3.2a,b) is even simpler than that used for the sea-surface temperature analysis. That is, the present model is a one layer rather than a two-layer transmission-emission model.

Equations (3.2a,b) as written show the two measured quantities I₃ and I₄ as depending on six unknowns: I_(G3), I_(G4), ε₃, ε₄, and T_(c). The values of I_(G3) and I_(G4) however, can be determined from nearby image pixels where there are no cirrus clouds present; and if ε₃ (v) and ε₄ (v) are replaced by band-averaged constant emissivities ε₃ and ε₄, there is a known functional relationship between the two emissivities as evidenced in the article by S. C. Ou, K. N. Liou, W. M. Gooch, and Y. Takano, "Remote Sensing of Cirrus Cloud Parameters Using Advanced Very-High-Resolution Radiometer 3.7- and 10-9 μm Channels", Applied Optics, Vol. 32, No. 12, 1993, pp. 2171-2180, incorporated herein by reference.

    (1-ε.sub.3)=(1-.sub.4).sup.m(T.sbsp.c.sup.)        (3.3a)

where ##EQU17## Function L(T_(c)) in equation (3.3b) is derived equation (3.1a) above. Now the forward model in equation (3.2a,b) can be approximated by

    I.sub.3 (v)=(1-ε.sub.4).sup.m(T.sbsp.c.sup.) ·I.sub.G3 (V)+[1-(1-ε.sub.4).sup.m(T.sbsp.c.sup.) ]·B(v,T.sub.c)(3.4a)

    I.sub.4 (v)=(1-ε.sub.4 ·I.sub.G4 (v)+ε.sub.4 ·B(v,T.sub.c)                                    (3.4b)

Such approximations are validated within the following articles incorporated herein by reference. Robert P. d'Entremont, Michael K. Griffin, and James T. Bunting, "Retrieval of Cirrus Radiative Properties and Altitudes Using Multichannel Infrared Data", Fifth Conference on Satellite Meteorology and Oceanography, Boston, Mass., American Meteorology Society, Sept. 3-7, 1990, is pp. 4-9; Robert P. D'ENTREMONT, Donald P. Wylie, J. William Snow, Michael K. Griffin, and James T. Bunting, "Retrieval of Cirrus Radiative and Spatial Properties Using Independent Satellite Data Analysis Techniques:, Proceedings Sixth Conference on Satellite Meteorology and Oceanography, Atlanta, Ga., American Meteorology Society, Jan. 5-10, pp. 17-20; Robert P. d'Entremont, Donald P. Wylie, Daniel C. Peduzzi, and Joseph Doherty, "Retrieval of Cirrus Radiative and Spatial Properties Using Coincident AVHRR and HIRS Satellite Data", Passive Infrared Remote Sensing of Clouds and the Atmosphere. SPIE-The International Society for Optical Engineering, Vol, 1934, 1993, pp. 180-196; S. C. Ou, K. N. Liou, W. M. Gooch, and Y. Takano, "Remote Sensing of Cirrus Cloud Parameters Using Advanced Very-High-Resolution Radiometer 3.7- and 10.9-μm Channels", Applied Optics, Vol. 32, No. 12, 1993, pp. 2171-2180; Robert P. d'Entremont and Gary B. Bustafson, "Support of Environmental Requirements for Cloud Analysis and Archive: Detection and Analysis of Cirrus Clouds Using Passive Infrared Satellite Data", Memo from Remote Sensing Group, Atmospheric and Environmental Research, Inc., 840 Memorial Drive, Cambridge, Mass., March 1996; and Algorithm Summaries, Robert d'Entremont, Jan. 22, 1997, VIIRS-40.4.6, pp. 35-37.

When equations (3.4a,b) are integrated over two bands, one in the 3.7 mm region and one in the 10-12 μm region, the four measured values of ∫I₃,dv, ∫I₄,dv, ∫I_(G3) dv, and ∫I_(G4) dv give two equations in two unknowns, providing enough information to use the statistical methodology of Section 1. For the 3.7 μm band, which we call band 3, equation (3.4a) is integrated from v_(3a) to v_(3b), and for the 10-μm band, which we call band 4, equation (3.4b) is integrated from v_(4a) to v_(4b). ##EQU18##

We next define two unknowns θ_(e) and θ_(c) using

    ε.sub.4 =ε.sub.04 +θ.sub.c           (3.6a)

    T.sub.c =T.sub.0 +θc                                 (3.6b)

where ε₀₄ and T₀ are typical (or average) values for, respectively, the cloud emissivity in band 4 and the cloud temperature. They are chosen so that the true values of the cloud temperature and emissivity lie reasonably close to T₀ and ε₀₄, making θ_(c) and θ_(e) small quantities.

The integrals of B(v,T_(c)) between v_(ia) and v_(ib) (for I=3,4) are then linearized with respect to θ_(c). ##EQU19## Equations (3.7a) come from a Taylor series expansion of B(v,T_(c)) about T_(c) =T₀, and equations (2.2a,b) are used to write the integrals in equations (3.7b,c). The next expression which has to be linearized is f(ε₄,T_(c)). We calculate

    f.sub.00 =f(ε.sub.04, T.sub.0)                     (3.8a) ##EQU20## with function L(T.sub.c) defined in equation (3.1a) above. Now f(ε.sub.4, T.sub.c) can be approximated as

    f(ε.sub.4,T.sub.c)≅f.sub.00 +f.sub.ε0 θ.sub.e +f.sub.TO θ.sub.c                     (3.8e)

using θ_(c) and θ_(e) as defined in equations (3.6a,b).

The last step in the linearization process requires the definition of two more variables θ_(G3) and θ_(G4).

    I.sub.G3 =I.sub.0G3 +θ.sub.G3                        (3.9a)

    I.sub.G4 =I.sub.0G4 +θ.sub.G4                        (3.9b)

In effect, the determination of I_(G3) and I_(G4) from nearby cloudless pixels is treated as two more measurements. The values of I_(0G3) and I_(0G4) are chosen to make θ_(G3) and θ_(G4) small, so they must be typical (or average) values of the I_(G3), I_(G4) radiances.

Equations (3.5a,b) and (3.9a,b) are combined to obtain the complete set of forward equations using (3.6b), (3.7a), and (3.8e) to represent ε₄, f, and the integrals over B.

    I.sub.3 =(f.sub.00 f.sub.ε0 θ.sub.e +f.sub.TO θ.sub.e)(I.sub.0G3 +θ.sub.G3)+[1-f.sub.00 -f.sub.ε0 θ.sub.e -f.sub.TO θ.sub.c ](b.sub.3 +β.sub.3 θ.sub.c)                                            (3.10a)

    I.sub.4 =(1-ε.sub.04 -θ.sub.e)(I.sub.0G4 +θ.sub.G4)+(ε.sub.04 +θ.sub.e)(b.sub.4 +β.sub.4 θ.sub.c)                                            (3.10b)

    I.sub.G3 =I.sub.0G3 +θ.sub.G3                        (3.10c)

    I.sub.G4 =I.sub.0G4 +θ.sub.G4                        (3.10d)

Equations (3.10a-d) are linearized by keeping terms of 0(θ_(c)), 0(θ_(e)), 0(θ_(G3)), 0(θ_(G4)); terms of 0(θ_(e) θ_(c)), 0(θ_(e),θ_(G3)), 0(θ_(e) θ_(G4)), 0(θ_(c) θ_(G3)), 0(θ_(c) θ_(G4)), 0(θ_(G3) θ_(G4)) are neglected, as are terms of 0(θ_(c) ^(n)), 0(θ_(e) ^(n)), 0(θ_(G3) ^(n)), 0(θ_(G4) ^(n)) for n>1. The method followed in the linearization process is that terms which are constant or linear in the θ's are kept, and terms which depend on higher powers of the θ's are dropped.

    I.sub.3 =θ.sub.c [ε(I.sub.0G3 -b.sub.3)+β.sub.3 (1-f.sub.00)]+θ.sub.e [f.sub.ε0 (I.sub.0G3 -b.sub.3)]+θ.sub.G3 [f.sub.00 ]+[f.sub.00 (I.sub.0G3 -b.sub.3)+b.sub.3 ]                                       (3.11a)

    I.sub.4 =θ.sub.c [β.sub.4 ε.sub.04 ]+θ.sub.e [b.sub.4 -I.sub.0G4 ]+θ.sub.G4 [1-ε.sub.04 ]+[I.sub.0G4 (1-ε.sub.04)+b.sub.4 ε.sub.04 ]           (3.11b)

    I.sub.G3 =I.sub.0G3 θ.sub.G3                         (3.11c)

    I.sub.G4 =I.sub.0G4 θ.sub.G4                         (3.11d)

Equations (3.11a-d) can be written in matrix form as ##EQU21## To check the accuracy of equations (3.12a-e) we represent bands 3 and 4 by the VIIRS point design for the 3.7 μm band and one of the 10-12 μm bands, as shown in Table 5. A noise vector w=(w₁,W₂,W₃,W₄)^(T) of zero-mean gaussian random numbers is added to the right-hand side of (3.12a).

    I=Hθ+s+w                                             (3.13a)

The noise amplitudes in bands 3 and 4 are σ₃ and σ₄ respectively (see Table 5 FIG. 13), so the covariance matrix for the w vector is ##EQU22## Equation (3.13b) assumes the I₃, I_(G3) measurements have the same σ₃ noise amplitude; it assumes that the I₄, I_(G4) measurements have the same σ4 noise amplitude; and it assumes the errors in all four measurements I₃, I₄, I_(G3), I_(G4) are uncorrelated. Equation (3.13a) has the form given in equation (1.8) above, so equation (1.9a) can be used to calculate C=(H^(T) C⁻¹ H)⁻¹ to get the expected error in our estimate of θ_(c), the cirrus cloud temperature. As will be shown below, and algorithm (see section 7 may be implemented which uses (1.9a) with the constants and matrices specified in Table 5 and equations (3.12e), (3.13b) to calculate .increment.T_(c) =0.836 deg K when T₀ =230 deg K, ε₀₄ =0.4, and the values of I_(0G3), I_(0G4) are the band 3 and 4 radiances of a black body at 290 deg K. The corresponding error in the ice-particle size of the cirrus cloud is given by ##EQU23## Table 6 (FIG. 14) is derived from the last page of Section 7 and shows that for a range of T_(c) values running between 220 deg K and 240 deg K the ice-particle size L runs between 41 μm and 88 μm, the absolute error in L is between 1.4 μm and 2.8 μm, and the relative error is between 2.9% and 3.3%. These values are a reasonable match to the threshold requirements of the ice-particle size EDR (5% relative error and 2 μm absolute error for ice particles between 0 and 50 μm in size, see reference 15). Again, as in Section 2, the calculations have produced believable results.

This statistical methodology can also be used to examine the effect of modeling error. If the forward model for predicting the values of I₃ and I₄ in equations (3.10a,b) does not predict exactly the I₃, I₄ atmospheric radiances, this can be regarded as another source of random error and include the amount by which the model is incorrect in the noise vector w. The modeling error can be treated as measurement "noise" because what is desired of the instruments to measure is the I₃, I₄ values predicted by the inaccurate model of (3.10a,b). Instead, what is obtained is the radiance coming from the unknown but accurate model, i.e. physical reality. There is no way of knowing whether the modeling error in any given measurement acts to compensate for the measurement noise, or to make the discrepancy worse, so σ_(mod) 3 and σ_(mod) 4 is taken as the modeling error in bands 3 and 4 respectively, to be independent sources of error obeying a zero-mean gaussian distribution. The covariance matrix C is now written as ##EQU24## The I_(G3), I_(G4) errors have not been changed because these are the radiances reaching the bottom of the cirrus cloud layer and are measured directly from nearby cloudless pixels. They cannot be affected by modeling error. They might, however, have somewhat different values in the nearby pixels from what they in fact are in the cirrus cloud pixels. We can, therefore, specify that σ_(G3) and σ_(G4) are the errors in I_(G3), I_(G4) respectively coming from the use of nearby cloudless pixels to estimate the values of I_(G3), I_(G4) in the cirrus cloud pixels. Including them in the formula for the covariance matrix C gives ##EQU25## Estimating reasonable values for a σ_(mod) 3, σ_(mod) 4, σ_(G3), σ_(G4) and then using (3.15) to specify matrix C in equation (1.9a) gives us a quick and easy way of including the modeling and procedure errors in the error budget for the ice-particle and cirrus temperature EDR's. Many of the other VIIRS-based EDR's can have their modeling and procedure errors including in their C matrices in a similar way.

Note also how easy it is to find a performance improvement of an EDR algorithm caused by the addition of an extra measurement to the forward model. Suppose the suggestion is made to add a 6.7 μm measurement, which shall be called band 5, to the cirrus EDR algorithm. Note that the relationship between ε₅, the cirrus emissivity in band 5, and ε₄, the cirrus emissivity in band 4 is given.

    (1-ε.sub.5)=(1-ε.sub.4).sup.n(T.sub.c)     (3.16)

where n(T_(c)) is a specified function of the cloud temperature T_(c). We define, following the procedure used on f(ε₄,T_(c)), that (see equations (3.8a-e))

    g(ε.sub.4,T.sub.c)=(1-ε.sub.4).sup.m(T.sbsp.c.sup.)(3.17a)

    g.sub.00 =g(ε.sub.04, T.sub.0)                     (3.17b) ##EQU26## Following the pattern of equation (3.11a), the two equations shown below are combined, equations (3.18a,b), with equations (3.11a-d) to get the new forward model.

    I.sub.5 =θ.sub.c [g.sub.TO (I.sub.0G5 -b.sub.5)+β.sub.5 (1-g.sub.00)]+θ.sub.e [g.sub.ε0 (I.sub.0G5 -b.sub.5)]+θ.sub.G5 [g.sub.00 ]+[g.sub.00 (I.sub.0G5 -b.sub.5)+b.sub.5 ]                                       (3.18a)

    I.sub.G5 =I.sub.0G5 θ.sub.G5                         (3.18b)

We next define b₅, β₅ in terms of the wavenumber limits v_(5a), v_(5b) of band 5. ##EQU27## The in-band signals for the band 5 radiances coming from pixels with and without the cirrus clouds are, respectively, I₅ and I_(G5). The value of I_(0G5) is chosen to make θ_(G5) small. Augmenting H, θ, and s to include equations (3.18a,b) gives ##EQU28## The matrix C corresponding to equation (3.13b) is ##EQU29## where σ₅ is the noise amplitude of the band 5 signal. Modeling and procedure errors can be added to (3.19d) the same way they were added to equation (3.15). Equation (1.9a) can now be used to find .increment.t_(c), the expected error in the cirrus cloud temperature, and equation (3.13c) then gives .increment.L, the expected error in the size of the cirrus ice particles.

Section 4: Computer-Based Forward Models

Some forward models are too complex to be easily written down algebraically. The models connecting atmospheric temperature and water-vapor profiles with the CrIs infrared spectra are one obvious example, and the VIIRS-based ocean-color and aerosol EDR's may be another (see reference 16). Typically, what does exist in these cases is a complex, computer-based forward model going from the EDR inputs to a collection of output radiances. Fortunately, this is all the information needed to find matrix H. We outline below how this is down, defining the variables in terms of the CrIS-based EDR's because the best-known forward models for these EDR's, such as MODTRAN or FASCODE, are entirely computerized.

The computer-based forward model is represented by the vector ##EQU30## where C is the total number of radiance outputs of the forward model. For the CrIS-based EDR's, it is the total number of FTIR channels in the measured infrared spectra (C is approximately 1500 for the current point design). The forward model inputs are represented by the vector ##EQU31## where S is the total number of these inputs. For the CrIS-based EDR's it is about twice the total number of atmospheric layers because each layer is characterized by both a temperature and a water-vapor concentration (S is typically about 70 for this type of forward model). An exact representation of the forward model can be written as

    f=f(T)                                                     (4.3a)

or

    f.sub.i =f.sub.i (T)=f.sub.i (T.sub.1, T.sub.2, . . . ,T.sub.5)(4.3b)

for I=1,2, . . . ,C.

The linearized forward equations are ##EQU32## where T₀ is the collection of input values for some standard atmosphere which is a good approximation for the true temperature and water-vapor profiles. ##EQU33## We calculate f₀, the vector of forward model outputs for the standard atmosphere T₀.

    f.sub.0 =f(T.sub.0)                                        (4.5a)

or ##EQU34## We define a new vector δT.sup.[k] to be an S-dimensional column vector which is zero everywhere except for the d'th component which is δT_(k). ##EQU35## We run the forward model for T₀ +δT.sup.[k], giving us a vector of outputs which we call f.sup.[k].

    f.sup.[k] =f(T.sub.0 +σT.sup.[k])                    (4.7a)

or ##EQU36## Vector f₀ represents the output of the computer-based forward model for the standard atmosphere inputs T₀, and f.sup.[k] is the output of the forward model with all the standard atmosphere inputs left unchanged except for the k'th input, which has been changed by an amount δT_(k). To create a complete collection of f.sup.[k] vectors the forward model must be run S times.

Returning to equation (4.4a) above, we recognize that it can be written as

    f=Hθ+f.sub.0                                         (4.8a)

where

    θ=T-T.sub.0,                                         (4.8b)

and H_(ij), the ij'th component of matrix H, is given by ##EQU37## Equations (4.8a-c) represent a restatement of equations (1.3b-d) in Section 1 above. The forward model predictions f in equation (1.8a) have no errors; but the actual CrIS spectra, of course, will. The measurement errors are represented in the usual way as an additive random vector w on the right-hand side of equation (4.8a), giving us

    f.sub.measured =Hθ+s+w                               (4.8d)

where s=f₀. The random vector w is characterized by a covariance matrix C coming from the point design of the CrIS. Equation (4.8d) is the same as equation (1.8), so as soon as matrix H is known, equation (1.9a) can be used to calculate the expected errors in the retrieved values of T₀ +θ. From (4.8b) we see that ##EQU38## where (f.sup.[k] -f₀)_(I) is the i'th component of the vector difference (f.sup.[k] -f₀). Equation (4.9) shows that the k'th column of H can be calculated by subtracting f₀, the standard atmosphere outputs of the forward model, from f[k], the forward-model outputs when only the k'th input of the standard atmosphere has been changed by δT_(k). This procedure requires the forward model to be run S+q times, once to calculate f₀ and S times to calculate f.sup.[k] for k=1,2, . . . ,S.

For the CrIS-based EDR's, such as the temperature and water-vapor profiles, there usually exists a matrix C.sub.θ describing the covariance of the expected EDR values between different levels of the atmosphere. Therefore, equation (1.11b) should be used to estimate the errors in these EDR's. In general, a covariance matrix C.sub.θ describes what is known about the statistics of the desired EDR's before any measurements are made. By making conservative or aggressive assumptions about C.sub.θ, we can determine conservative or aggressive values for the expected errors in these Cr1S-based EDR's. When no assumptions are made about the EDR values, the C.sub.θ matrix can be set to zero and equation (1.9a) used to find the expected errors in the retrieved EDR's. Equation (4.9) can be used to find H for any computer-based forward model; accordingly, the above detailed steps are also applicable to VIIRS-based EDR's.

Section 5

Suppose we have two straight lines written in functional form as f₁ (v) and f₂ (v) where

    f.sub.1 (v)=m.sub.1 v+b.sub.1,

    f.sub.2 (v)=m.sub.2 v+b.sub.2,

and we want to see whether

    f.sub.1 (v)≅A·f.sub.2 (v)

for some constant A. One quick way of doing this is to see whether ##EQU39## because

    m.sub.1 v+b.sub.1 =A·(m.sub.2 v+b.sub.2)

can only be true (for all values of v) when m₁ =Am₂ and b₁ =Ab₂, or ##EQU40##

Section 6

This algorithm is used to calculate the expected error in the sea surface temperature measured by satellite IR radiometers operating in the 800 to 1000 cm⁻¹ atmospheric spectral window and utilizes formulas on pages 97 and 391 of Steven M. Kay's textbook Fundamentals of Statistical Signal Processing: Estimation Theory, PTR Prentice Hall, Englewood Cliffs, N.J. 1993, ISBN 0-13-345711-7.

Start by defining constants used in the calculations. Variable n is the power of v (wavenumber) over temperature used to match (exp(x)-1) (-1) in the Planck function for wavenumbers between 800 and 100 cm⁻¹ and temperatures between 310 and 290 deg K. Variable CB is a constant used to normalize the measured inband radiances when applying the model equations. It has units of watt-cm2/sr/(cm-deg K) n. The atmospheric transmittance is modeled as 1-K(v)*A where K(v) is a straight line of slope M in units of cm. (This transmittance model assumes that the radiometer measurements are made over a wide enough band that we can neglect the rapid up-and-down variation of the atmospheric transmittance with respect to v.)

N:=4.4

CB:=1.98·10⁻¹²

M:=8.0·10⁻⁴

Define the F4 and F5 functions which will be used to calculate the elements of matrix H as defined on pages 97 and 391 of Steven Kay's text. Variables va and vb are, respectively, the beginning and ending wavenumbers of the spectral band being measured by the satellite radiometer. ##EQU41## Specify two bands for the radiometer measurement, band 1 between 800 and 870 cm⁻¹ and band 2 between 885 and 971 cm⁻¹. Variable σ1 is the standard deviation of the measurement error in band 1 in units of watt/cm2/sr, and σ2 is the standard deviation of the measurement error in band 2 in units of watt/sm2/sr. We set up matrices C and H as defined on pages 97 and 391 of Steven Kay's book. ##EQU42## Utilize the formula on page 97 of Steven Kay's book to calculate the expected error in the estimate of TS n, the n'th power of the sea-surface temperature in units of degK n. Matrix Cu is the covariance matrix for the estimate of Ts n and the unknown atmospheric transmittance parameter. ##EQU43## The corresponding error in the measured sea-surface temperature is σTs, its standard deviation in degK. It is expected that the sea-surface temperature lies between 290 and 310 degK, so we calculate a table of corresponding σTs values.

    ______________________________________                                         Ts.sub.0 : = 290. Ts.sub.1 : = 295. Ts.sub.2 : = 300. Ts.sub.3 : = 305.        Ts.sub.4 : = 310. i: = 0..4                                                     ##STR1##                                                                       ##STR2##                                                                              Ts.sub.i                                                                            σTs.sub.i                                                   ______________________________________                                                 290  0.7672964                                                                 295  0.7239716                                                                 300  0.6837606                                                                 305  0.6463934                                                                 310  0.6116269                                                         ______________________________________                                    

Section 7

This algorithm is used to estimate the error involved in finding the ice-particle size of cirrus clouds using VIIRS measurements (based on the 1997 NPOES proposal) in the 3.7 micron (band 3) and 10.9 micron (band 4) radiance bands. The analysis presented below is based on information about the dependence of the bands 3 and 4 radiances in a paper "Remote Sensing of Cirrus Cloud Parameters Using Advanced Very-High-Resolution Radiometer 3.7 and 10.9 μm Channels" in Applied Optics, Vol. 32, No. 12, April 1993 by S. C. Ou, K. N. Liou, W. M. Gooch, and Y. Takanco, pp. 2171-2180, (see especially eqns. (1a,b) on p. 2172 and eqn. (4) on p. 2173), and incorporated herein by reference. The statistical analysis is based on eqns. (4.30) and (4.32) on page 97 of Fundamentals of Statistical Signal Processing: Estimation Theory by Steven M. Kay, PTR Prentice Hall, Englewood Cliffs, N.J. 1993 (ISBN 0-13-345711-7).

Create a function called PLnck(x) which calculates 1(exp(x)-1). When x is close to zero it replaces exp(x)-1 and returns 1/x exactly. ##EQU44## PLnck(x):=if(|x|>0.01, PLnck(x), PLnck2(x)) Generate C1 and C2 black body constants used to specify Planck's curve. Variable C1 has units of watt/sm2/sr/(cm-1)4 and variable C2 has units of degK/cm-1.

    C1:=1.191062·10.sup.12 C2:=1·438786

Specify Ibb₋₋ Winvcminvsr as a function of wavenumber a (in cm-1) and temperature T (in degK). The units of Ibb₋₋ Winvminvsr are Watts/cm/sr. ##EQU45## Specify IdT₋₋ WinvcminvsrinvdegK, the derivative of Ibb₋₋ Winvcminvsr, as a function of wavenumber σ (in cm-1) and temperature T (in degK). The units of IdT₋₋ WinvcminvsrinvdegK are W/cm/sr/degK. ##EQU46## Specify the beginning and ending wavenumbers (in cm-1) of band 3, called v3a and v3b respectively, and the beginning and ending wavenumbers of band 4, called v4a and v4b respectively.

v3a₋₋ invcm:=2545

v3b₋₋ invcm:=2817

v4a₋₋ invcm:=885

v4b₋₋ invcm:=971

Variable nedt3₋₋ degK is the NEDT (in degK) expected for the band 3 point design of VIIRS in the 1997 NPOES proposal. Variable nedt4₋₋ degK is the NEDT (in degk) for the band 4 point design of VIIRS in the 1997 NPOES proposal. Tnedt3₋₋ degK and Tnedt4₋₋ degK are the band 3 and band 4 scene temperatures at which the NEDT values are evaluated.

Nedt3₋₋ degK:=0.2

Tnedt3₋₋ degK:=271.0

nedt4₋₋ degK:=0.1

Tnedt4₋₋ degK:=293.0

Variable σ3₋₋ Winvcm2invsr is the expected measurement error (standard deviation of the additive noise) in band 3 in units of Watt/sm2/sr, and σ4₋₋ Winvcm2invsr is the expected measurement error (standard deviation) in band 4 in units of Watt/sm2/sr. ##EQU47## Calculate the constants a3₋₋ Winvcm2invsr, β3₋₋ Winvcm2invsrinvdegK and a4₋₋ Winvcm2invsr, β4₋₋ Wincm2invsrinvdegK which are used to linearize with respect to temperature the integrals over bands 3 and 4 of the Planck black-body function. The integrals are linearized about temperature TO₋₋ degK (in units of degK). TO₋₋ degK:=230.0 ##EQU48## Define function L₋₋ μm m, the cirrus cloud ice-particle size (in microns) as a function of cloud temperature Tc (in deg K). The best fit constants c0₋₋ μm, c1₋₋ μminvdegK, c2₋₋ minvdegK2, and c3₋₋ minvdegK3 come from curve fitting of cloud data. C0₋₋ μm:=326.3 c1₋₋ μminvdegK:=12.42 c2₋₋ μminvdegK2:=0.197 c3 μminvdegK3:=0.0012 L₋₋ μm(tc):=c0μm+c1₋₋ μminvdegK·(Tc-273)+c2₋₋ μminvdegK2 ·(Tc-273)² +c3₋₋ μminvdegK3 (Tc-273)³

Define function m (dimensionless), the exponent used in the forward radiance model for bands 3 and 4. It depends on L₋₋ μm, the cirrus cloud ice-particle size in microns. The constants a0s, a2s₋₋ μm, and A2s₋₋ μm2 define the dependence of m on L₋₋ μm. a0s:=0.722 a1s₋₋μm:= 55.08 a2s₋₋ μm2:=-174.12 ##EQU49## Define the derivative of function m with respect to its argument L₋₋ μm, mdL₋₋ invμm (in units of μm-1) ##EQU50## Define the derivative of L₋₋ μm with respect to its argument Tc, LdT₋₋ μminvdegK (in units of μm/degK) LdT₋₋ μminvdegK(Tc):=c1₋₋ μminvdegK+2·c2₋₋ μminvdegK2·(Tc-273)+3·c3₋₋ μminvdegK3·(Tc-273)²

Define the derivative of m with respect to the cirrus cloud temperature Tc, mdT₋₋ invdegK (in degK-1).

MdT₋₋ invdegK(Tc):=mdL₋₋ invμm(L₋₋ μm(Tc))·LdT₋₋ μminvdegK(Tc)

Specify ε04, the dimensionless average emissivity of the cirrus clouds. ε04:=0.4

Define the values of the constants f00, f.sub.ε, 0, and fT0₋₋ invdegK.

F00:=(1-ε04)^(m) (T0₋₋ degK)

F00=0.586773

fε0:=-m(T0₋₋ degK)·(1-ε04)^(m)(T0.sbsp.--^(degK))-1

fε0=1.020631

fT0₋₋ invdegK:=(1-ε04)^(m)(T0.sbsp.--^(degK)) ·In(1-ε044)·mdT₋₋ invdegK(T0₋₋ degK)

fT0₋₋ invdegK=0.003404

Create typical values for the amount of in-band radiance in bands 2 and 3 reaching the cirrus clouds from the ground by estimating the surface temperature as Ts₋₋ degK (in degK) and integrating the Planck function at this temperature over bands 3 and 4. Variable IOG3₋₋ Winvcm2invsr and IOG4₋₋ Winvcm2invsr are these in-band radiances for bands 3 and 4 respectively (in units of watts/cm2/sr). ##EQU51## Define the H matrix used on page 97 of Steven Kay's text (referenced at start of this spreadsheet).

H11:=fT0₋₋ invdegK-(I0G3₋₋ Winvcm2invsr-a3₋₋ Winvcm2invsr)+β3₋₋ Winvcm2invsrinvdegK·(1-f00)

H12:=fε0·(I0G3₋₋ Winvcm2invsr-a3₋₋ Winvcm2invsr)

H13:=f00

H21:=β4₋₋ Winvcm2invsrinvdegK·ε04

H22:=a4₋₋ Winvcm2invsr-I0G4₋₋ Winvcm2invsr

H24:=I-ε04 ##EQU52## Specify the noise covariance matrix C specified in page 97 of Steven Kay's text (referenced at the beginning of the spreadsheet) ##EQU53## Calculate the covariance matrix of the errors in the determination of unknown parameters: the cirrus cloud temperatures, the emissivity of the cirrus in band 4, the deviation of the radiance reaching the cloud from the ground from its nominal values in band 3, and the same radiance deviation in band 4. Here we are only interest in the error in the estimate of the cirrus cloud temperature, whose variance will be given by the "1,1" value of the error covariance matrix Cθ.

    Cθ:=(H.sup.T ·C.sup.-1 ·H).sup.-1 ##EQU54## The standard deviation of the error in the cirrus temperature (in degK) is σTc.sub.-- degK. ##EQU55## Use the previously-defined function L.sub.-- μm, which depends on the cirrus cloud temperature in degK, to create a table of the error in L.sub.-- μm expected given the above error σ Tc.sub.-- degK in the estimated cloud temperatures.

    ______________________________________                                         i: = 0,,21 Tc.sub.i : = 220 + (i - 1) · 1 Lmicron.sub.i : =           L.sub.-- μm(Tc.sub.i)                                                       Lmicron.sub.-- err.sub.i : = LdT.sub.-- μminvdegK(Tcp.sub.i)                · σTc.sub.-- degK                                                ##STR3##                                                                      Lmicron.sub.i Lmicron.sub.-- err.sub.i                                                                  Rel.sub.-- errL.sub.i                                 ______________________________________                                         41.1152       1.373158   0.033398                                              42.7606       1.380519   0.032285                                              44.4184       1.393903   0.031381                                              46.0958       1.413309   0.03066                                               47.8          1.438738   0.030099                                              49.5382       1.47019    0.029678                                              51.3176       1.507664   0.029379                                              53.1454       1.55116    0.029187                                              55.0288       1.60068    0.029088                                              56.975        1.656222   0.029069                                              58.9912       1.717786   0.029119                                              61.0846       1.785374   0.029228                                              63.2624       1.858984   0.029385                                              65.5318       1.938616   0.029583                                              67.9          2.024271   0.029813                                              70.3742       2.115949   0.030067                                              72.9616       2.213649   0.03034                                               75.6694       2.317372   0.030625                                              78.5048       2.427118   0.030917                                              81.475        2.542886   0.031211                                              84.5872       2.664677   0.031502                                              87.8486       2.79249    0.031788                                              ______________________________________                                    

Conclusion

The statistical methodology outlined in Section 1, and applied to an algebraic forward model for the sea-surface temperature in Section 2 and an algebraic forward model for the size of cirrus ice particles in Section 3 provides a measure of the accuracy in which a satellite may yield parameters of the earth's surface. In Section 2 the predicted EDR performance matched NOAA's past experience for two-band sea-surface temperature retrieval, and in Section 3 the predicted performance matched NOAA's threshold requirements. Thus, in both cases, the statistical methodology led to reasonable results. The performance of an augmented EDR algorithm can also be determined when an extra band is used to provide more spectral information, as well as how the effect of modeling error could be easily included in the covariance matrix describing the measurement uncertainties. Section 4 showed how to obtain the matrix elements of the linearized forward model (which the statistical methodology requires) when the forward model is only available as a computer program. The forward model for the CrIS-based EDR's, as well as some of the more complicated VIIRS-based EDR's, may in fact, only exist in computerized form.

Current techniques for deriving VIIRS-based EDR's seem to rely on creating forward models consisting of N equations in N unknowns and then solving for one or more of the unknowns to create the desired algorithm. Clearly, such a procedure cannot be applied to forward models consisting of N equations in M unknowns when N>M; yet the information provided by the extra equations might well lead to more accurate EDR algorithms. This sort of situation, where the forward model has more equations than unknowns, is exactly the sort of situation which is easily handled by the statistical methodology according to the present invention. The current VIIRS-based EDR algorithms also have no way of using correlations between the M unknowns of the forward mode, and this is additional extra information which can be easily handled by this methodology. While EDR algorithms based on the methodology will have come from linearized forward models and so may be less accurate than EDR algorithms based on the full non-linear forward models, their ability to use extra information may provide a performance advantage over their non-linear counterparts, especially when the desired EDR's exhibit only moderate variations about a climatic mean. For these reasons, algorithms derived from equations (1.9b) and (1.11 c,d) in Section 1 may well provide a significant advance over the current state of the art in VIIRS-based EDR retrieval.

Moreover, the statistical method according to the present invention provides a quick and easy means of estimating the expected performance of an EDR algorithm. Current methods for determining performance envisions a series of Monte-Carlo, end-to-end simulations based on combinations of EDR forward models with the simulated sensors and the proposed EDR algorithms. The inputs of the forward models would be changed randomly and compared again and again to the outputs coming from the EDR algorithms. This is a truly ambitious program of number crunching and may, in the end, provide only imprecise knowledge of the EDR algorithms' performance, especially if only a relatively few Monte-Carlo runs are performed. Although some type of end-to-end Monte-Carlo simulation will be required to verify performance predictions achieved by other means (such as the statistical methodology discussed in Sections 1-4 above), this sort of simulation is at best a clumsy way of examining different algorithms for the same EDR. As we have shown in Sections 1-4, the present statistical methodology only needs an EDR's linearized forward model to predict the EDR algorithm's behavior; compared to an end-to-end simulation it well more quickly and easily suggest which version of an EDR algorithm has the best performance. When it becomes time to perform an end-to-end simulation, this methodology will, by having already examined the algorithm's behavior, provide an approximate answer to the simulation's results, helping to avoid gross errors in the computer calculations. 

What is claimed is:
 1. A method of estimating expected errors of environmental data parameters based on radiance measurements obtained from visible infra red radiometric satellite sensors (VIIRS) orbiting the earth, comprising the steps of:obtaining N radiance measurements of a surface body I_(I), . . . I_(N) defining a matrix I depending on p unknown surface and atmospheric parameters T_(i), . . . T_(p), defining a matrix T; generating a forward model I=f(T) for obtaining said I radiance measurements from said p parameters, where f(T) comprises an N element matrix; choosing an initial set of values for said p parameters and linearizing f(T) about said initial values to obtain a linearized forward model I=s+Hθ as f_(i) (T₀₁, T₀₂, . . . T_(0p))+ΣH_(ij) θ_(j) where I=1,2, . . . N and θ_(j) =T_(j) -T_(oj) and H_(ij) =δ_(i) f/δT_(j) and where θ is a column matrix ₁ θ, . . . _(p) θ and H is a matrix of H_(ij) values; adding measurement noise vector w of noise values to said forward model; determining the covariance of said measurement noise w to obtain a covariance matrix C; and manipulating the matrices H and C according to the equation C_(EDR=)(H^(T) C⁻¹ H)⁻¹ to obtain a matrix element of C_(EDR) indicative of the expected errors in said values of T₁. . . T_(p) parameters.
 2. The method according to claim 1, wherein said measurement noise w comprises N random numbers w₁, . . . ,w_(N).
 3. The method according to claim 2, wherein said N random numbers have a jointly gaussian distribution probability.
 4. The method according to claim 3, wherein said jointly gaussian distribution probability is a zero mean gaussian distribution.
 5. The method according to claim 1, wherein said matrix element C_(EDR) specify expected errors in the values θ₁, . . . θ_(p).
 6. The method according to claim 5, wherein the expected errors in the values θ₁, θ₂, . . . θ_(p) represent the differences in said p surface and atmospheric parameters from said initial set of values.
 7. A method of predicting errors in sea surface temperature measurements based on radiance measurements I obtained from visible infra red radiometric satellite sensor (VIIRS) orbiting the earth having at least a first bandwidth and a second bandwidth corresponding to a particular window region, comprising the steps of:applying a forward model for determining sea-surface temperature from said satellite radiance measurements I over said particular wavelength window region as I(v,atm)=B(v,T_(s))·t(V,atm)+B(v,T_(a))·[1-t(v,atm)] where T_(s) is measured sea-surface temperature, Ta represents a temperature indicative of atmospheric self-radiance in said wavelength window region, t represents a parameter indicative of transmissivity from said sea surface to the top of the atmosphere, v is a wavenumber in said window region, and atm is representative of atmospheric parameters determining the value of said transmittance at wavenumber v, and where B(V,T) is a Planck function; approximating said transmissivity parameter t(v,atm) as 1-A_(atm) ·K(v) in said window region where K is a function only of wavenumber and where A_(atm) is indicative of atmospheric parameters; performing a straight line approximation of K(v) to obtain K_(L) (V); evaluating B(V,T) over a limited range of V and T parameter values to obtain I(V,atm)=C_(B) is a constant and θ_(s) =(T_(s))^(n), and θ_(a) =A_(atm) (T_(a) ^(n) -T_(s) ^(n)); measuring noise amplitudes over said first and second bandwidths to obtain a first noise amplitude σ₁, and a second noise amplitude σ₂ ; integrating said radiance function measurements I(v,atm) over said first and second bandwidths to obtain said radiance I as a function of a matrix H and a matrix θ such that I=Hθ, where I is a matrix function of radiances I1 and I2 corresponding to said first and second bandwidths, and θ is a matrix representation of θ_(s) and θ_(a), and where H is a function of V; adding to said radiance I a noise matrix w corresponding to said noise amplitudes to obtain a radiance I=Hθ+w; determing the covariance of said noise matrix w to obtain a covariance matrix C;and manipulating the matrices H and C according to the equation (H^(T) C⁻¹ H)⁻¹ to obtain a matrix C_(EDR) having an element indicative of the expected variance in said sea surface temperature variable T_(s) ^(n).
 8. The method according to claim 7, wherein said wavelength window region is 10-12 μm.
 9. The method according to claim 8, wherein said T_(a) varies by less than 1 Kelvin within said 10-12 μm window.
 10. The method according to claim 9, wherein said limited range over which B(V,T) is evaluated further comprises a temperature range of between 285 K and 315 K.
 11. The method according to claim 10, wherein said limited range over which B(V,T) is evaluated further comprises a wavenumber range of between 800 cm⁻¹ and 1000 cm⁻¹.
 12. The method according to claim 11, wherein said Planck function B(V,T) is C₁ V³ P(C₂ VT⁻¹) where C₁ is 1.191×10⁻¹² watt/cm² /sr/(cm⁻¹)⁴, and where P(C₂ VT⁻¹)=1/(e^(c2VT-1) -1).
 13. The method according to claim 7, wherein the step of approximating t(v,atm) as 1-A_(atm) ·K(v) is applied only to wide bandwidth radiance signals.
 14. A method estimating an error involved in determining ice particle size of cirrus clouds using radiance measurements obtained from visible infrared radiometric satellite sensors (VIIRS) in a first and second radiance bands, said VIIRS sensors orbiting the earth and obtaining said radiance measurements from a body located on a surface of said earth, comprising the steps of:creating a Planck function P(x) such that the value of said P(x) is 1/x when x approaches 0, and is 1/(e^(x) -1) otherwise; determining first and second black body constants C1 and C2 associated with specifying a curve of said Planck function P(x); determining a first parameter as a function of wavenumber σ and temperature τ; determining a second parameter as a function of wavenumber σ and temperature τ; determining beginning and ending wavenumbers associated with said respective first and second radiance bands; storing an expected temperature error NEDT expected for said first band associated with the design of said VIIRS sensor and a second expected temperature error NEDT4 expected for said second band associated with the design of said VIIRS sensor; storing first and second temperature measurement values TNEDT and TNEDT4, respectively, with each of said first and second respective radiation bands over which said expected temperature errors NEDT and NEDT4 are evaluated; obtaining an expected measurement error σ₃ in said first band by integrating said first parameter as a function of said first expected temperature error and first temperature measurement value over said first radiance band and obtaining a second expected measurement error σ₄ associated with said second radiance band by integrating said first parameter as a function of said second expected temperature error and said second temperature measurement value over said second radiance band; determining a set of constants A3, B3, A4, and B4 for linearizing with respect to temperature integrals over said first and second bands of said Planck black body function P(x), wherein said constant A3 comprises integrating said first parameter as a function of temperature T₀ over said first radiance band; wherein said constant B3 is evaluated by integrating said second parameter at the temperature T₀ over said first radiance band; wherein said A4 constant evaluated by integrating said first parameter at said temperature T₀ over said second radiance band; and wherein said constant B4 is evaluated by integrating said second parameter at said temperature T₀ over said second radiance band; defining a cirrus cloud ice particle size parameter L as a function of cloud temperature T_(c) based on a plurality of constants c0, c1, c2 and c3 derived from a curve fitting of cloud data; defining a function M as a function of said L cirrus cloud ice particle size and constant values A0, A1, and A2 for associating said dependence of M on L; defining a derivative dM of function M with respect to L and a derivative dL of L with respect to temperature T_(c) ; defining a derivative dTM of M with respect to cirrus cloud temperature T_(c) ; determining an average emissivity of said cirrus clouds as ε₀₄ ; defining values associated with constants f₀₀, f.sub.ε0, and fT0 where f₀₀ is characterized by the function f₀₀ =(1-ε₀₄)^(m)(To)-1 and f.sub.ε0 is characterized by the function -m(T₀) (1-ε₀₄)^(m)(To)-1 and fT0 is characterized by the function FT0=(1-ε₀₄)^(m) To ·1n(1ε₀₄)·dTM(T₀); estimating the surface temperature T_(s) and integrating said Planck function P(x) at said surface temperature Ts over each of said first and second radiance bands to obtain a first estimated radiance I_(G3) associated with said first band, and a second estimated radiance I_(G4) associated with said second band; defining a matrix H having corresponding matrix elements indicative of said values fT0, f.sub.ε0, and f₀₀ ; determining a noise covariance matrix C based on said first and second expected measurement errors σ₃ and σ₄ ; determining a covariance matrix of errors in the estimation of said cirrus cloud temperature C.sub.θ where C.sub.θ is characterized by the equation C.sub.θ =(H^(T) ·C⁻¹ ·H)⁻¹ ; determining a standard deviation σT of said estimated errors in Cθ; and applying said σT to said L(Tc) particle size to obtain said estimate of error in said particle size.
 15. The method according to claim 14, further comprising the step of creating a table of errors L given said above error estimate σT in said estimated cloud temperatures using said previously defined function L which depends on said cirrus cloud temperature.
 16. The method according to claim 14, wherein said temperature measurements are measured in Kelvins.
 17. The method according to claim 14, wherein said first band comprises a 3.7 μm region and where said second band comprises a 10.9 μm region radiance band.
 18. The method according to claim 17, wherein said beginning and ending wave numbers in said first band are 2545 cm⁻¹ and 2817 cm⁻¹ respectively.
 19. The method according to claim 18, wherein said beginning and ending wavenumbers for said second band comprise 885 cm⁻¹ to 971 cm⁻¹ respectively.
 20. The method according to claim 19, wherein said first and second temperature measurement values are respectively 271 K and 293 K. 